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Abstract 

We apply Maxwell and Cattaneo's relaxation approaches to the analysis of strong Shockwaves 
in a two-dimensional viscous heat-conducting fluid. Good agreement results for reasonable values 
of Maxwell's relaxation times. Instability results if the viscous relaxation time is too large. These 
relaxation terms have negligible effects on slower-paced subsonic problems, as is shown here for 
two-roll and four-roll Rayleigh-Benard flow. 
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I. INTRODUCTION 



In 1867 James Clerk Maxwell^ noted that an initial shear stress in a dilute gas, (like 
air) when unsupported by an underlying shear motion, will decay with a relaxation time 
r = [f]/ P) (about 200 picoseconds for air), where r] is the shear viscosity and P the pressure. 
His governing relaxation equation for the shear stress modifies Newton's a = rje to read 

a + Ta = Tje . 

Here a is the stress, t] the viscosity, and e the strain rate. The superior dots represent 
comoving time derivatives. 

Nearly a century later Carlo Cattanec-^ argued that Fourier's law for heat conduction 
should be similarly modified, in order to avoid the supersonic heat flow implied by a parabolic 
(diffusion equation) transport law. One could equally well argue that a heat flux, when 
unsupported by a temperature gradient, would decay with a microscopic relaxation time r 
like Maxwell's. Cattaneo's approach can be written in a form like Maxwell's, but with a 
partial (fixed in space) rather than a comoving time derivative: 

Q + T{dQ/dt) = -kVT . 

Cattaneo's rationale for using a partial time derivative rather than one fixed in the material 
is unclear. Here Q is the heat flux, T the temperature, and n the heat conductivity. With 
Cattaneo's relaxation assumption, "heat waves" can propagate at about the speed of sound^. 
On physical grounds Maxwell's approach, with the comoving time derivative, seems more 
"realistic" than Cattaneo's. Cattaneo's form for the relaxation time makes no contribution 
at all in stationary steady-state problems such as the structure of a steady fluid shockwave. 

Oddly enough, modern treatments of time delay^i^ often use Cattaneo's partial-derivative 
formulation rather than Maxwell's comoving time derivative. The purpose of the present 
work is to elucidate the usefulness of the relaxation concept and to explore its limits in 
applications of fluid mechanics. In the following Sections we consider the relatively fast- 
paced steady shockwave problem as well as the slower-paced steady convective Rayleigh- 
Benard flow. A final Section summarizes our findings. For simplicity we use units in which 
the Boltzmann constant and atomic mass are both equal to one. 
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II. STRONG DENSE-FLUID SHOCKWAVES 



The structure of strong Shockwaves has long served as a testing ground for continuum 
models like the Navier- Stokes- Fourier equations (here given for a two-dimensional fluid with 
vanishing bulk viscosity, rjy = 0): 

p = -pV -v ; pv = -V ■ P ; pe = -Vv : P-V -Q ; 

P = /[Pcq + r/V ■ v] - r][Vv + Vv^] ; Q = -kVT . 

The time derivatives, here as before indicated by the superior dot, are all comoving deriva- 
tives, like Maxwell's, time rates of change in a coordinate frame moving with the fluid 
velocity v. Solving the three differential equations for the density p, velocity v, and energy e 
requires a knowledge of the pressure tensor P and heat flux vector Q. The simplest models 
are shown here, with two transport coefficients, the Newtonian shear viscosity rj and the 
Fourier heat conductivity k defined in the usual way. / is the unit tensor, with I^^ = lyy = 1 
and Ixy = lyx — 0. 

Landau and Lifshitz' analytic solution of the shockwave structure for a gas with con- 
stant transport coefficients and a shockwidth A provides a useful initial condition for both 
macroscopic continuum and microscopic molecular dynamics simulations^: 

P{^) = ^_x/x ^ e+x/A ^ { v{x),Pxx{x),Qx{x) } . 

Their solution smoothly interpolates the density between cold fluid, with density pc, and 
hot fluid, with p^. 

Molecular dynamics shockwave simulations^- have been carried out in the two different 
ways shown in Figure 1: (1) by following the two moving waves generated by the inelastic 
collision of two blocks of material; (2) by studying the single stationary wave formed with 
two boundary "treadmills" - on the left boundary cold fluid is introduced at the "shock 
speed" Vs while at the right boundary hot fluid is extracted at the slower speed Vg — Vp, 
where Vp is the "particle speed". In either case, in a coordinate frame centered on the 
shockwave the mass, momentum, and energy fluxes are all constant: 

{pv, Pxx + pv^, pv[e + {Pxx/ p) + (i'^/2)] + Qx} constant for all x . 
For "weak" shocks the Navier- Stokes- Fourier description is "good"-^^. For "stronger" shocks 
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(twofold compression) several contradictions to this simple description arise^-^^. To illustrate 
these points typical mechanical and thermal shockwave profiles are shown in Figure 2. 

First, the local longitudinal and transverse temperatures differ, often by more than a fac- 
tor of two (see Figure 2). Second, as is also shown in Figure 2, the shear stress {Pyy — Pxx)/'^ 
and the heat flux both lag behind the velocity gradient [dv^/dx) and the tempera- 
ture gradients {dT^x/dx) and (dTyy/dx), suggesting the presence of Maxwell- type relaxation 
time a ^ — . Third, the fact that temperature is so very anisotropic makes it necessary to 
consider separate xx and yy contributions to the heat fiux-^-i^: 



Fourth, the same anisotropicity also suggests including asymmetric divisions of the work 
and heat contributions (indicated by d) to the thermal energy change: 



Here Cy is the heat capacity per unit mass. Fifth, a mechanism for the decay of temperature 
anisotropy must also be included: 



Last, the molecular dynamics results imply that a bulk viscosity rjv, approximately equal 
to the shear viscosity, must be includedi^.. Though a continuum model incorporating all of 
these ideas is necessarily relatively complex, a successful implementation of all six of these 
additions to the Navier-Stokes-Fourier model is described in References^i,— , and^. 

In those works all of the continuum field variables were derived from molecular dynamics 
simulations using a short-ranged repulsive pair potential. 



The prefactor (lO/vr) was chosen to give a potential energy integral of unity for a random 
particle distribution at unit density: 




{pCv/2) T,, D [-aVv : Pxhermal - /3V ■ Q] ; 



{pCv/2) Tyy D [-(1 - a)Vv : Pxhermal - (1 - /3)V ■ Q] . 




0(r < 1) = (10/7r)(l -r)^ . 
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The initial zero-pressure zero-temperature state was compressed twofold to obtain a hot 
dense fluid state. Lucy's normalized weighting functioni^ii^ was used to compute spatial 
averages of the various fleld variables: 

rh 

w{r <h) = (5/7r/i^)[l - {r/h)f\l + 3(r//i)] ^ / 2Tirw{r) = 1 . 

Jo 

The smooth-particle average of the particle quantity fj is given by a weighted sum, 

(p(r)/(r)) = ^mjfjw{r - rj) ; p(r) = ^mjU;(r - rj) . 
j j 

This smooth-particle deflnition has two advantages: (1) all of the fleld variables deflned in 
this way have two continuous space derivatives; (2) the continuity equation (with fj equal 
to the particle velocity vj) is satisfled exactly: 

{ p = ^ rrijwir — rj) ; pt> = ^ mjVjw{r — rj) } — > p = — pV ■ v . 
j j 

Here p and pv are deflned everywhere in this way, not just at the particle locations. The 
range h of the "weighting function" w{r < h) is typically chosen so that about 20 particles 
contribute to fleld-point averages. With this approach the microscopic pressure tensor and 
heat flux vector at any point in space are expressed in terms of nearby individual particle 
contributions to these nonequilibrium fluxesi^i^S. 

To appreciate the effect of the various modiflcations of the Navier- Stokes- Fourier model we 
next study the stability of solutions using a continuum model which is a rough representative 
of the molecular dynamics results^"—. 

III. STABILITY STUDIES WITH AN IDEALIZED GRUNEISEN MODEL 

For stability studies we choose an equilibrium equation of state based on Griineisen's 
separation of the energy and pressure into cold and thermal parts: 

Peq = pe= (pV2) + 2pT ; e = (p/2) + 2T . 

A Shockwave satisfying all the conservation laws results when a cold fluid is compressed 
to twice its initial density by a shockwave moving toward that fluid at twice the particle 
velocity {vg = 2vp = 2). In this case the constant mass, momentum, and energy fluxes are 
respectively 

{pv = 2; P^^ + pv^ = (9/2) ; pv[e + {P^^/p) + (t;V2)] + = 6 } . 
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The various hydrodynamic variables then cover the following ranges within the Shockwave: 



[ 2 > v{x) >!];[!< p{x) < 2 ] ; [ (1/2) < e{x) < (5/4) ] ; 

[ (1/2) < Peq < (5/2) ] ; [ < < (1/8) ] . 

(Note that T^x can exceed the "hot" value of (1/8) within the Shockwave.) The details 
of the Shockwave structure depend upon the nonequilibrium constitutive relations for the 
shear stress and the heat flux. Next, we summarize two separate situations, (1) vanishing 
conductivity with a scalar temperature; (2) tensor conductivity, with separate longitudinal 
and transverse temperatures, with different contributions from work and heat. Both these 
models lead to the conclusion that the mechanical relaxation time cannot be too large. By 
contrast, the thermal relaxation time can be either "small" or "large" . 

A. Relaxation Without Heat Conduction 

The simplest case results when both heat conductivity and thermal anisotropy are omit- 
ted. Then the density and energy can both be eliminated from the three flux equations, 

pv = 2; (pe) -a + 2v = (9/2) ; 2[e + e - {v/2)a + (^72)] = 6 , 

giving the shear stress, 

^ — {Pyy Pxx) /'^ ~ Pxx i 

as a function of velocity: 

(7 = (3/'f;)(v-l)(v-2) < . 

Evidently the viscous stress is everywhere negative (compressive). If we introduce Maxwell's 
idea of comoving stress relaxation, 

(7 -\- Ta = (7 -\- Tv{da/dx) = a + Tv{da/dv){dv/dx) = r){dv/dx) , 

we find that the velocity gradient (dv/dx) diverges unless the ratio (r/r/) is sufficiently small: 

Ta < (ry/3) . 

It is physically reasonable that too long a memory can lead to instability in fast-paced 
complex flows like Shockwaves. On the other hand the relaxation equation by itself, with a 
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smooth strain increment localized near zero time {t = 0), 

1 

a + Ta = - — — - , 

[e-* + e+*] 

provides smooth solutions even for large t—^. The present analytic shockwave limit on 
To- < (^/3) is in full accord with two kinds of numerical simulations. First, the stationary 
flux equations can be solved for the temperature and stress fields, just as was indicated 
above for the case of vanishing conductivity. Second, it is possible to solve the dynamical 
equations for 

{ (dp/dt), (dv/dt), (de/dt) } or { (dp/dx), (dv/dx), (de/dx) } 

starting with the Landau-Lifshitz profile. The two methods agree. They show that the 
stress relaxation time in Shockwaves must be sufficiently small, t„ < {i]/3) for stability. 

We next extend the thermal constitutive model to include tensor temperature with 
anisotropic heat conduction. We also include separate relaxation times for the longitu- 
dinal and transverse heat fluxes, and the separation of work and heat into longitudinal and 
transverse parts^»^. 



B. Lack of Relaxation Without Viscosity 

Viscosity, as opposed to heat conduction, is essential to the shock process. To appreciate 
this need, consider the conservation equations for our simple model without viscosity and 
with the heat conductivity equal to unity: 

e = (p/2) + 2T ; pv = 2 ; pe + pv"^ = (9/2) ; pv[2e + {v^/2)] - (dT/dx) = 6 . 

According to the first three equations the temperature has its maximum value of (T^ax = 
0.238 > Thot = 0.125) within the shock: 

{p, V, T} = {1.4436, 1.38545, 0.23800} for T = T^ax ■ 

But the fourth (energy-fiux) equation gives (dT/dx) = 0.7106 for that thermodynamic state, 
contradicting the presence of a maximum. Thus this model, lacking viscosity, cannot sustain 
a stationary shockwave. 
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Exactly this same conclusion follows also for the inviscid ideal gas, with twofold compres- 
sion from unit density, pressure, and temperature, with Vg — and the wholly thermal 
pressure P = pe = pT. Because heat conductivity in the absence of viscosity is not enough to 
provide a shockwave, the relaxation effects are quite different for conductivity and viscosity, 
as we show next. 

C. Relaxation with Tensor Temperature, Apportioned Work and Heat 

The analysis becomes more complicated when heat flow is included, along with relax- 
ation and separated contributions of the heat and work to the longitudinal and transverse 
temperatures. Here the heat flux evolves following the tensor relaxation equation: 

Qx ~l~ '^qQx — ^Xx{dTxx/ Kyy{dTyy / (Ix) . 

The divergence of the heat flux provides net heating and is apportioned between the longi- 
tudinal and transverse temperatures: 

pT^x ^ -PidQJdx) ; pfyy D (1 - p){dQjdx) . 

The contributions of the heat flux divergence V • Q to heating are indicated by the inclusion 
symbol, "d". We include also an analogous separation of the thermodynamic work into 
longitudinal and transverse parts: 

pTxx ^ -a-Prhermal : Vv ; pfyy D (1 - Q;)PThermal : Vf . 

Finally, the two temperatures necessarily relax toward one another: 

TxX ^ i'^yy Txx) /TQ , Tyy D (Txx ^yy) / Tq . 

For simplicity we choose the two thermal relaxation times [ for the heat flux Q and the tem- 
perature anisotropicity {T^x — Tyy) ] to have a common value, tq. For illustrative purposes 
we emphasize the difference between the two temperatures by choosing the apportionment 
parameters a and /3 both equal to unity, so that both the work and the heat provide longi- 
tudinal heating, with the transverse temperature lagging behind. 

Then straightforward (at least for a computer) algebra provides solutions of the shockwave 
problem and reveals not one, but two restrictions on tq. For stable solutions to exist we 
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found in this way that the thermal relaxation time must be either sufficiently small or 
sufficiently large. Setting the distance scale of the shockwave with the constant transport 
coefficients 

T] 2,Kxx '^^yy ' 
computer algebra gives the following restrictions on the relaxation times: 

< r, < (1/3) ; TQ < (1/8) or tq > (1/4) . 

Figures 3 and 4 show typical continuum profiles using these constitutive relations. The 
continuum profiles were generated in two quite different ways: (1) solving the time-dependent 
equations for {p, v, e, a, Q} starting with the Landau-Lifshitz approximation; (2) solving the 
stationary flow equations for the mass, momentum, and energy fluxes using a computer 
algebra program (we used "Maple"). The latter approach provides page- long formulae for 
(du/dx), (dTrcx/dx), and (dTyy/dx) as well as numerical solutions. The stationary equations 
for the shockwave profile have no solution if the relaxation time for the shear stress is 
greater than (rj/S) or if the relaxation time for the heat flux lies between (k/8) and (k/4). 

To summarize, our findings for Shockwaves establish that momentum-fiux relaxation has 
to be "fast" for stability. Thermal relaxation can either be likewise fast or quite slow, with 
a window of instability separating these two regimes. Where the thermal relaxation is slow 
the shockwave structure is dominated by viscosity rather than conductivity. 

It is natural to speculate on the effect of relaxation in ordinary hydrodynamic situations. 
In order to see what consequences arise from these effects in subsonic fluid mechanics we 
next introduce delay into the hydrodynamic equations describing a compressible, conducting, 
viscous flow, the Rayleigh-Benard problem. 

IV. RAYLEIGH-BENARD FLOW 

To investigate the stability of moderate flows to the presence of viscous and thermal 
relaxation we revisit some finite-difference Rayleigh-Benard simulations of two-roll, four- 
roll, and six-roll flows2ii2^. The simulations picture a viscous conducting fluid, heated from 
below in the presence of a vertical gravitational field. Sufficiently strong heating causes a 
transition from static heat conduction to one of a number of nonequilibrium steady states 
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with stationary convection rolls. Stationary and transient sample flows are shown in Figures 
5 and 6. 

For the Rayleigh-Benard model we study here (equal kinematic viscosity and thermal 
diffusivity) the transition from static Fourier conduction to two-roll convection occurs near 
a Rayleigh number R of 1750: 

R = g{d\nV/dT)pH^AT/{i^DT) = H'^/{vDt) . 

The fluid is confined to a rectangular box, periodic on the sides, with the gravitational con- 
stant g = (l/H) chosen to give constant density in the nonconvecting case. H is the height 
of the cell, equal here to half the width. AT is the difference between the hot temperature 
at the base {T^ = 1.5) and the cold temperature at the top of the cell {Tc = 0.5). u 
and Dt (chosen equal, for convenience) are the kinematic viscosity and thermal diffusivity 
(both with units of [length"^ /time]) . For simplicity we choose all values of the relaxation 
times equal and do not distinguish between the longitudinal and transverse temperatures, 
Txx = Tyy. Our model continuum fluid obeys the ideal gas equation of state: 

Peq = pT = pe ] 7]v = ] r] = 2Kxx = '^l^yy = 1 ■ 

Numerical results for this model are given as a function of Rayleigh number in 
References^ andr^. Simulations with the various relaxation times all equal to 0.1 repro- 
duced this earlier work perfectly. As an example, the two-roll problem of References^ and^, 
with a Rayleigh number of 40,000 gives per-cell kinetic energies of (K^/N) + (Ky/N) = 
0.00373 + 0.00357. We carried out many special cases with a Rayleigh Number of 40,000, 
which produces stationary steady states. Whether two-roll or four-roll solutions are obtained 
is sensitive to the initial conditions^. We began with a very weak two-roll velocity field as 
the initial condition in an if x ly box with the coordinate origin at its center: 

Vx (X sm{2'7ix /W) sm{2ny / H) ; Vy oc cos(27ra;/iy) cos^ny/H) . 

We found solutions for t = brj = 5k, and r = lOrj = 10k but instability when r was doubled 
again to 20t] = 20k. These additions of relaxation to the Navier- Stokes- Fourier equations 
lowered the horizontal kinetic energy and raised the vertical, with both effects on the order 
of parts per thousand. Thus relaxation in subsonic flows has only relatively small effects in 
the regime of stable solutions. 
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V. CONCLUSIONS 



Molecular dynamics simulations have established the facts that delay times on the order 
of a collision time, as envisioned by Maxwell, affect shockwave structure in a substantial way. 
Cattaneo's approach, with partial time derivatives, has no effect on shockwave structure. 
Shockwaves are dominated by viscosity, so that stress relaxation must be relatively rapid. 
Thermal relaxation, important for chemical relaxation, can be either fast or slow. 

In ordinary subsonic fluid mechanics the effects of time delays are relatively small. As 
a result, thermal anisotropicity is ordinarily ignored in continuum mechanics. It is a sub- 
stantial effect in shocks, with repercussions for chemical reaction rates. In our continuum 
simulations we have assumed relaxation equations with comoving time derivatives, 

a + Tr^d = rit , Q + tqQ = -kVT , 

rather than partial derivatives. If a were replaced with (da/dt) there would be no relaxation 
at all in a stationary problem like the shockwave and Rayleigh-Benard problems studied here. 

The Maxwellian relaxation times cause no trouble solving conventional moderate flow 
problems like Rayleigh-Benard convection. The problem areas suggested by this work in- 
clude (1) formulating optimum choices for locally- averaged hydrodynamic variables with 
the general goal of maximizing the accuracy of macroscopic descriptions of microscopic re- 
sults and (2) developing theoretical models for the estimation of the relaxation parameters 
measured in the dynamical simulations. 

A logical approach to problem (1) above would use "entropy production" as a tool^^. 
In the Rayleigh-Benard problem entropy production is proportional to the squares of the 
nonequilibrium fluxes, o"^ and Q^. If these are computed locally, with a weight function 
w{r < h) then h can be chosen such that the internal entropy production matches the 
boundary sources and sinks of entropy. Evidently too small/large an h gives too large/small 
an entropy production, so that h can be chosen to be "just right". Problem (2) would have 
to begin with some nonequilibrium simulations tailored to the direct measurement of delay 
and relaxation. 

Finally, the presence of delay has some pedagogical importance. Delay in the results 
of time-reversible motion equations (molecular dynamics) breaks the time-symmetry which 
would otherwise lead to a logical contradiction between time-reversible molecular dynamics 
and conventional irreversible continuum mechanics^. 
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I. FIGURE CAPTIONS 

• Figure 1. Two colliding fluid blocks generate symmetric Shockwaves (velocities ±fp) 
as the blocks, moving at ±[vs — Vp] collide and come to a stop (shown above). Two 
treadmill boundaries, one fast (velocity Vg) and one slow (velocity Vs — Vp), maintain 
a single stationary shockwave in the center of the system (shown below). 

• Figure 2. Density, pressure, internal energy, temperature tensor, and heat flux in a 
strong shockwave simulation using molecular dynamics. In the cold unshocked material 
the nearest-neighbor spacing is unity. The hot shocked fluid has a density exactly twice 
that of the cold material. Figure based on data described in reference^^. 

• Figure 3. Solution of the continuum model for twofold compression with the Griineisen 
equation of state using = (1/10) and tq = tr = 1. The mass, momentum, and en- 
ergy fluxes are {2, (9/2), 6} Compare with Figure 4 noting particularly the differences 
between T^x and Tyy. 

• Figure 4. Solution of the continuum model for twofold compression with the Griineisen 
equation of state using = tq = r/j = (1/10). The mass, momentum, and energy 
fluxes are {2, (9/2), 6} Compare with Figure 3. 

• Figure 5. Transient flow fleld for the Rayleigh-Benard problem at time = 1000. The 
initial state was two weakly-rotating rolls. The viscous, heat-conducting, compressible 
fluid is heated at the bottom and cooled at the top with a gravitational field directed 
downward. The vertical boundaries at the sides are periodic. The number of compu- 
tational cells shown here is 80 x 40 = 3200. The transport coefficients, rj = k, = (1/5) 
were selected to give a Rayleigh number of 40,000. The relaxation times were set equal 
to unity: = tq = 1. 
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• Figure 6. A fully-converged four-roll structure evolved from the flow field shown in 
Figure 5. Here the time is 10,000. The Rayleigh number is 40,000 and the viscosity 
and heat conductivity, rj = k = (1/5), have equal relaxation times, = tq = 1. 
The final kinetic energy is (KJN) + (Ky/N) = 0.001144 + 0.004133 = 0.005277. The 
number of computational cells is A'^ = 80 x 40. 
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Figure 1 
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Figure 3 
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Figure 4 
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Figure 5 
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Figure 6 
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